/**
 * TODO:
 * @Implement Stentiford Thinning Method
 * @Implement Morphological Thinning
 * @Implement 
 * 
 */

#include <util/common_macros.h>

#include "../../matrix/matrix.h"
#include "immorph.h"


// <editor-fold defaultstate="collapsed" desc="Auxiliar functions to Zhang and Suen Method">

inline static int __checkN8(real * data, int cols) {
	int n8 = 0;

	if (*data == 0) return 0;

	// data[-1][-1]
	data = data - cols - 1;
	if (*(data++) != 0) n8++;
	if (*(data++) != 0) n8++;
	if (*(data) != 0) n8++;

	// data[0][-1]
	data += cols - 2;
	if (*(data++) != 0) n8++;
	data++; // skip central pixel
	if (*(data) != 0) n8++;

	// data[1][-1]
	data += cols - 2;
	if (*(data++) != 0) n8++;
	if (*(data++) != 0) n8++;
	if (*(data) != 0) n8++;

	return n8;
}

inline static int __checkNT(real * data, int cols) {
	int nt = 0;

	data = data - cols;
	if (*data == 0 && *(data + 1) > 0) nt++;
	data = data + 1;
	if (*data == 0 && *(data + cols) > 0) nt++;
	data = data + cols;
	if (*data == 0 && *(data + cols) > 0) nt++;
	data = data + cols;
	if (*data == 0 && *(data - 1) > 0) nt++;
	data = data - 1;
	if (*data == 0 && *(data - 1) > 0) nt++;
	data = data - 1;
	if (*data == 0 && *(data - cols) > 0) nt++;
	data = data - cols;
	if (*data == 0 && *(data - cols) > 0) nt++;
	data = data - cols;
	if (*data == 0 && *(data + 1) > 0) nt++;

	return nt;
}// </editor-fold>

/**
 * This function implements the thinning algorithm of Zhang and Suen described
 * at article:
 * 
 * T. Y. Zhang and C. Y. Suen. 1984. A fast parallel algorithm for thinning
 * digital patterns. Commun. ACM 27, 3 (March 1984), 236-239.
 * DOI=10.1145/357994.358023 http://doi.acm.org/10.1145/357994.358023
 * 
 * @param rows - number of rows
 * @param columns - number of columns
 * @param mat - input matrix pointer
 * @param thin - output matrix pointer
 * @return number of erodes.
 */
int thinning_zhang_suen(int rows, int columns, real * mat, real * thin) {

	int r, c, _r, i;
	int size = rows * columns;
	int np, sp;
	int count1 = 0;
	int count2 = 0;
	real p2, p4, p6, p8;

	real tmp[size];

	// <editor-fold defaultstate="collapsed" desc="Step 01">
	for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
		for (c = 1; c < columns - 1; c++) {

			i = _r + c;

			// = skip zero pixels
			if (mat[i] == 0) {
				tmp[i] = 0;
				continue;
			}
			// = get connectivity N(p) of p
			np = __checkN8(&mat[i], columns);

			// = Checks the First Condition (2 <= N(p) <= 6) ===================
			if (2 <= np && np <= 6) {

				// = get number of transitions 0->1 around p
				sp = __checkNT(&mat[i], columns);

				// = Checks the Second Condition (S(p) == 1) ===================
				if (sp == 1) {
					p2 = *(&mat[i] - columns);
					p4 = *(&mat[i] + 1);
					p6 = *(&mat[i] + columns);

					// = Checks the Third Condition (p2 * p4 * p6 == 0) ========					
					if (p2 == 0 || p4 == 0 || p6 == 0) {
						p8 = *(&mat[i] - 1);

						// = Checks the Fourth Condition (p4 * p6 * p8 == 0) ===
						if (p4 == 0 || p6 == 0 || p8 == 0) {

							tmp[i] = 0;
							count1++;
							continue;
						}
					}
				}
			}
			tmp[i] = 1;
		}
	}
	// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Step 02">
	for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
		for (c = 1; c < columns - 1; c++) {

			i = _r + c;

			// = skip zero pixels ==============================================
			if (tmp[i] == 0) {
				thin[i] = 0;
				continue;
			}
			// = get connectivity N(p) of p ====================================
			np = __checkN8(&tmp[i], columns);

			// = Checks the First Condition (2 <= N(p) <= 6) ===================
			if (2 <= np && np <= 6) {
				sp = __checkNT(&tmp[i], columns);

				// = Checks the Second Condition (S(p) == 1) ===================
				if (sp == 1) {
					p2 = *(&tmp[i] - columns);
					p4 = *(&tmp[i] + 1);
					p8 = *(&tmp[i] - 1);

					// = Checks the Third Condition (p2 * p4 * p8 == 0) ========					
					if (p2 == 0 || p4 == 0 || p8 == 0) {
						p6 = *(&tmp[i] + columns);

						// = Checks the Fourth Condition (p2 * p6 * p8 == 0) ===
						if (p2 == 0 || p6 == 0 || p8 == 0) {

							thin[i] = 0;
							count2++;
							continue;
						}
					}
				}
			}
			thin[i] = 1;
		}
	}
	// </editor-fold>

	return (count1 + count2);
}

/**
 * This function implements the thinning algorithm of Lu and Wang described
 * at article:
 * 
 * H.E. Lu and P.S.P. Wang. 1985. An improved fast parallel algorithm for
 * thinning digital patterns. Proc. of the IEEE Conf. on computer vision and
 * pattern recognition. pp. 364-367.
 * 
 * This method is a little modification of Zhang and Suen method.
 * 
 * @param rows - number of rows
 * @param columns - number of columns
 * @param mat - input matrix pointer
 * @param thin - output matrix pointer
 * @return number of erodes.
 */
int thinning_lu_wang(int rows, int columns, real * mat, real * thin) {

	int r, c, _r, i;
	int size = rows * columns;
	int np, sp;
	int count1 = 0;
	int count2 = 0;
	real p2, p4, p6, p8;

	real tmp[size];

	// <editor-fold defaultstate="collapsed" desc="Step 01">
	for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
		for (c = 1; c < columns - 1; c++) {

			i = _r + c;

			// = skip zero pixels
			if (mat[i] == 0) {
				tmp[i] = 0;
				continue;
			}
			// = get connectivity N(p) of p
			np = __checkN8(&mat[i], columns);

			// = Checks the First Condition (2 <= N(p) <= 6) ===================
			if (3 <= np && np <= 6) {

				// = get number of transitions 0->1 around p
				sp = __checkNT(&mat[i], columns);

				// = Checks the Second Condition (S(p) == 1) ===================
				if (sp == 1) {
					p2 = *(&mat[i] - columns);
					p4 = *(&mat[i] + 1);
					p6 = *(&mat[i] + columns);

					// = Checks the Third Condition (p2 * p4 * p6 == 0) ========					
					if (p2 == 0 || p4 == 0 || p6 == 0) {
						p8 = *(&mat[i] - 1);

						// = Checks the Fourth Condition (p4 * p6 * p8 == 0) ===
						if (p4 == 0 || p6 == 0 || p8 == 0) {

							tmp[i] = 0;
							count1++;
							continue;
						}
					}
				}
			}
			tmp[i] = 1;
		}
	}
	// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Step 02">
	for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
		for (c = 1; c < columns - 1; c++) {

			i = _r + c;

			// = skip zero pixels ==============================================
			if (tmp[i] == 0) {
				thin[i] = 0;
				continue;
			}
			// = get connectivity N(p) of p ====================================
			np = __checkN8(&tmp[i], columns);

			// = Checks the First Condition (2 <= N(p) <= 6) ===================
			if (2 <= np && np <= 6) {
				sp = __checkNT(&tmp[i], columns);

				// = Checks the Second Condition (S(p) == 1) ===================
				if (sp == 1) {
					p2 = *(&tmp[i] - columns);
					p4 = *(&tmp[i] + 1);
					p8 = *(&tmp[i] - 1);

					// = Checks the Third Condition (p2 * p4 * p8 == 0) ========					
					if (p2 == 0 || p4 == 0 || p8 == 0) {
						p6 = *(&tmp[i] + columns);

						// = Checks the Fourth Condition (p2 * p6 * p8 == 0) ===
						if (p2 == 0 || p6 == 0 || p8 == 0) {

							thin[i] = 0;
							count2++;
							continue;
						}
					}
				}
			}
			thin[i] = 1;
		}
	}
	// </editor-fold>

	return (count1 + count2);
}

/**
 * Image thinning / skeletonization
 * @param Mat - input matrix
 * @param Thin - output matrix. The skeletonized image
 * @param n - number of iterations. If n is less than 0 the function runs until there is no more to thinning.
 */
int imp_image_thinning_zhang_suen(ImpMatrix * Mat, ImpMatrix ** Thin, int n) {

	int it = 0;
	int rows = Mat->rows;
	int columns = Mat->columns;
	ImpMatrix * Tmp = imp_matrix_new_copy(Mat);
	int changes;

	if (*Thin == NULL) (*Thin) = imp_matrix_new(Mat->rows, Mat->columns);

	real * tmp = Tmp->data;
	real * thin = (*Thin)->data;

	changes = 1;

	while (n-- && changes) {
		changes = thinning_zhang_suen(rows, columns, tmp, thin);
		SWAP(tmp, thin, real *);
		it++;
	}

	imp_object_delete(Tmp);

	return it;
}

/**
 * Image thinning / skeletonization
 * @param Mat - input matrix
 * @param Thin - output matrix. The skeletonized image
 * @param n - number of iterations. If n is less than 0 the function runs until there is no more to thinning.
 */
int imp_image_thinning_lu_wang(ImpMatrix * Mat, ImpMatrix ** Thin, int n) {

	int it = 0;
	int rows = Mat->rows;
	int columns = Mat->columns;
	ImpMatrix * Tmp = imp_matrix_new_copy(Mat);
	int changes;

	if (*Thin == NULL) (*Thin) = imp_matrix_new(Mat->rows, Mat->columns);

	real * tmp = Tmp->data;
	real * thin = (*Thin)->data;

	changes = 1;

	while (n-- && changes) {
		changes = thinning_lu_wang(rows, columns, tmp, thin);
		SWAP(tmp, thin, real *);
		it++;
	}

	imp_object_delete(Tmp);

	return it;
}

int hit_kernel3x3(int rows, int columns, real * mat, short int kernel[9], int * hit) {

	int m;
	int r, c, _r, i;
	real p;
	int count = 0;

	for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
		for (c = 1; c < columns - 1; c++) {
			i = _r + c;

			if (mat[i] == 0) continue;

			boolean h = TRUE;

			for (m = 0; m < 9; m++) {
				if (kernel[m] != 0 && kernel[m] != 1) continue;

				if (m == 0) p = mat[i - columns - 1 ];
				else if (m == 1) p = mat[i - columns ];
				else if (m == 2) p = mat[i - columns + 1 ];
				else if (m == 3) p = mat[i - 1 ];
				else if (m == 4) p = mat[i];
				else if (m == 5) p = mat[i + 1 ];
				else if (m == 6) p = mat[i + columns - 1 ];
				else if (m == 7) p = mat[i + columns ];
				else if (m == 8) p = mat[i + columns + 1 ];

				if (p > 0) p = 1;

				if (p != kernel[m]) {
					h = FALSE;
					break;
				}
			}

			if (h) {
				hit[i] = 1;
				count++;
			}
		}
	}
	return count;
}

int thinning_stentiford(int rows, int columns, real * mat, real * thin) {

	short kernel[4][9] = {
		{
			2, 0, 2,
			2, 1, 2,
			2, 1, 2
		},
		{
			2, 2, 2,
			1, 1, 0,
			2, 2, 2
		},
		{
			2, 1, 2,
			2, 1, 2,
			2, 0, 2
		},
		{
			2, 2, 2,
			0, 1, 1,
			2, 2, 2
		}

	};

	int k;
	int i;
	int count = 0;
	int size = rows*columns;

	int * tmp = mem_allocn(size, int);

	for (i = 0; i < size; i++) {
		thin[i] = mat[i];
		tmp[i] = 0;
	}

	for (k = 0; k < 4; k++) {

		count += hit_kernel3x3(rows, columns, thin, kernel[k], tmp);

		for (i = 0; i < size; i++) {
			if (tmp[i]) {
				tmp[i] = 0;
				thin[i] = 0;
			}
		}
	}

	mem_free(tmp);

	return count;
}

int imp_image_thinning_stentiford(ImpMatrix * Mat, ImpMatrix ** Thin, int n) {

	int it = 0;
	int rows = Mat->rows;
	int columns = Mat->columns;
	ImpMatrix * Tmp = imp_matrix_new_copy(Mat);
	int changes;

	if (*Thin == NULL) (*Thin) = imp_matrix_new(Mat->rows, Mat->columns);

	real * tmp = Tmp->data;
	real * thin = (*Thin)->data;

	changes = 1;

	while (n-- && changes) {
		changes = thinning_stentiford(rows, columns, tmp, thin);
		SWAP(tmp, thin, real *);
		it++;
	}

	imp_object_delete(Tmp);

	return it;
}

int medial_axis_transform(int rows, int columns, real * mat, real * thin) {

	/*
		short kernel[8][9] = {
			{
				2, 0, 2,
				2, 1, 2,
				2, 1, 2
			},
			{
				2, 2, 0,
				2, 1, 2,
				1, 2, 2
			},
			{
				2, 2, 2,
				1, 1, 0,
				2, 2, 2
			},
			{
				1, 2, 2,
				2, 1, 2,
				2, 2, 0
			},
			{
				2, 1, 2,
				2, 1, 2,
				2, 0, 2
			},
			{
				2, 2, 1,
				2, 1, 2,
				0, 2, 2
			},
			{
				2, 2, 2,
				0, 1, 1,
				2, 2, 2
			},
			{
				0, 2, 2,
				2, 1, 2,
				2, 2, 1
			}
		};
	 */
	/*
		short kernel[4][9] = {
			{
				2, 0, 2,
				2, 1, 2,
				2, 1, 2
			},
			{
				2, 2, 2,
				1, 1, 0,
				2, 2, 2
			},
			{
				2, 1, 2,
				2, 1, 2,
				2, 0, 2
			},
			{
				2, 2, 2,
				0, 1, 1,
				2, 2, 2
			}
		
		};
	 */

	short kernel[8][9] = {
		{0, 0, 0, 2, 1, 2, 1, 1, 1},
		{2, 0, 0, 1, 1, 0, 2, 1, 2},
		{1, 2, 0, 1, 1, 0, 1, 2, 0},
		{2, 1, 2, 1, 1, 0, 2, 0, 0},
		{1, 1, 1, 2, 1, 2, 0, 0, 0},
		{2, 1, 2, 0, 1, 1, 0, 0, 2},
		{0, 2, 1, 0, 1, 1, 0, 2, 1},
		{0, 0, 2, 0, 1, 1, 2, 1, 2}
	};

	int k;
	int i;
	int count = 0;
	int size = rows*columns;

	int * tmp = mem_allocn(size, int);

	for (i = 0; i < size; i++) {
		thin[i] = mat[i];
		tmp[i] = 0;
	}

	for (k = 0; k < 8; k++) {

		count += hit_kernel3x3(rows, columns, thin, kernel[k], tmp);

		for (i = 0; i < size; i++) {
			if (tmp[i]) {
				tmp[i] = 0;
				thin[i] = 0;
			}
		}
	}

	mem_free(tmp);

	return count;
}

int imp_image_morph_medial_axis_transform(ImpMatrix * Mat, ImpMatrix ** Thin) {

	int it = 0;
	int rows = Mat->rows;
	int columns = Mat->columns;
	ImpMatrix * Tmp = imp_matrix_new_copy(Mat);
	int changes;

	if (*Thin == NULL) (*Thin) = imp_matrix_new(Mat->rows, Mat->columns);

	real * tmp = Tmp->data;
	real * thin = (*Thin)->data;

	changes = 1;

	while (changes) {
		changes = medial_axis_transform(rows, columns, tmp, thin);
		SWAP(tmp, thin, real *);
		it++;
	}

	imp_object_delete(Tmp);

	return it;
}

int thin_features(int rows, int columns, real * mat, real * thin) {

	short kernel[8][9] = {
		{
			2, 0, 2,
			2, 1, 2,
			2, 1, 2
		},
		{
			2, 2, 0,
			2, 1, 2,
			1, 2, 2
		},
		{
			2, 2, 2,
			1, 1, 0,
			2, 2, 2
		},
		{
			1, 2, 2,
			2, 1, 2,
			2, 2, 0
		},
		{
			2, 1, 2,
			2, 1, 2,
			2, 0, 2
		},
		{
			2, 2, 1,
			2, 1, 2,
			0, 2, 2
		},
		{
			2, 2, 2,
			0, 1, 1,
			2, 2, 2
		},
		{
			0, 2, 2,
			2, 1, 2,
			2, 2, 1
		}
	};
	
	int k;
	int i;
	int count = 0;
	int size = rows*columns;

	int * tmp = mem_allocn(size, int);

	for (i = 0; i < size; i++) {
		thin[i] = mat[i];
		tmp[i] = 0;
	}

	for (k = 0; k < 8; k++) {

		count += hit_kernel3x3(rows, columns, thin, kernel[k], tmp);

		for (i = 0; i < size; i++) {
			if (tmp[i]) {
				tmp[i] = 0;
				thin[i] = 0;
			}
		}
	}

	mem_free(tmp);

	return count;
}

int imp_image_thin_features(ImpMatrix * Mat, ImpMatrix ** Thin, int n) {

	int it = 0;
	int rows = Mat->rows;
	int columns = Mat->columns;
	ImpMatrix * Tmp = imp_matrix_new_copy(Mat);
	int changes;

	if (*Thin == NULL) (*Thin) = imp_matrix_new(Mat->rows, Mat->columns);

	real * tmp = Tmp->data;
	real * thin = (*Thin)->data;

	changes = 1;

	while (n-- && changes) {
		changes = thin_features(rows, columns, tmp, thin);
		SWAP(tmp, thin, real *);
		it++;

	}

	imp_object_delete(Tmp);

	return it;
}

int thin_kernel3x3(int rows, int columns, real * mat, short int kernel[9], real * thin) {

	int m;
	int r, c, _r, i;
	real p;
	int count = 0;

	for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
		for (c = 1; c < columns - 1; c++) {
			i = _r + c;

			if (mat[i] == 0) {
				thin[i] = 0;
				continue;
			}

			thin[i] = 1;

			boolean hit = TRUE;

			for (m = 0; m < 9; m++) {
				if (kernel[m] != 0 && kernel[m] != 1) continue;

				if (m == 0) p = mat[i - columns - 1 ];
				else if (m == 1) p = mat[i - columns ];
				else if (m == 2) p = mat[i - columns + 1 ];
				else if (m == 3) p = mat[i - 1 ];
				else if (m == 4) continue;
				else if (m == 5) p = mat[i + 1 ];
				else if (m == 6) p = mat[i + columns - 1 ];
				else if (m == 7) p = mat[i + columns ];
				else if (m == 8) p = mat[i + columns + 1 ];

				if (p > 0) p = 1;

				if (p != kernel[m]) {
					hit = FALSE;
					break;
				}
			}

			if (hit) {
				thin[i] = 0;
				count++;
			}
		}
	}
	return count;
}
